# Figures: uses data file created in stata
library(foreign)
library(readstata13)
library(ggplot2)
setwd("C:/Users/Tom/Dropbox (Personal)/Papers/Swiss Panel Data/")
a=read.dta13("SHP_finaldataR.dta")

#####################################################################

# Figure 1: Unemployment Risk

## 1. Age ##
pdf(file="fig1a.pdf")

d=c("18-24", "25-39", "40-64")  
b=c("Unemployment Rate","Mean subjective Unemployment Risk")  
corrs=matrix(nrow=2,ncol=3)
dimnames(corrs)=list(b,d)

corrs[1,1:3]=c(7.12,3.63,2.73)
corrs[2,1]=(mean(a$subj_urisk[which(a$age<25)]))-1
corrs[2,2]=(mean(a$subj_urisk[which(a$age<40&a$age>24)]))-1
corrs[2,3]=(mean(a$subj_urisk[which(a$age<65&a$age>39)]))-1                  

par(mfrow=c(1,2))
barplot(corrs[1,],col=c("gray60","gray60","gray60"),las=1,ylab="Unemployment Rate (Per Cent)",
        ylim=c(0,8),xlab="",cex.lab=1.4,cex.axis=1.5,cex.names=1.2)
barplot(corrs[2,],col=c("gray60","gray60","gray60"),
        ylab="Mean Subjective Unemployment Risk",ylim=c(0,1.2),yaxt="n"
        ,xlab="",cex.lab=1.4,cex.names=1.2)
axis(2,at=c(0,0.2,0.4,0.6,0.8,1,1.2),labels=c(1,1.2,1.4,1.6,1.8,2,2.2),las=1,cex.axis=1.5)
mtext("(a) Age Groups",side=3,outer=T,line=-2,cex=1.7)
dev.off()

####

## 2. Education ##
a$educat=as.numeric(a$educat)
pdf(file="fig1b.pdf")

d=c("Level 1", "Level 2", "Level 3")  
b=c("Unemployment Rate","Mean subjective Unemployment Risk")  
corrs=matrix(nrow=2,ncol=3)
dimnames(corrs)=list(b,d)

corrs[1,1:3]=c(6.61,3.48,2.38)
corrs[2,1]=(mean(a$subj_urisk[which(a$educat<11)]))-1
corrs[2,2]=(mean(a$subj_urisk[which(a$educat>11&a$educat<19)]))-1
corrs[2,3]=(mean(a$subj_urisk[which(a$educat==19)]))-1                  

par(mfrow=c(1,2))
barplot(corrs[1,],col=c("gray60","gray60","gray60"),las=1,
        ylab="Mean Unemployment Rate (Per Cent)",
        ylim=c(0,7),xlab="",cex.lab=1.4,cex.axis=1.5,cex.names=1.1)
barplot(corrs[2,],col=c("gray60","gray60","gray60"),
        ylab="Mean Subjective Unemployment Risk",ylim=c(0,1.2),yaxt="n"
        ,xlab="",cex.lab=1.4,cex.names=1.1)
axis(2,at=c(0,0.2,0.4,0.6,0.8,1,1.2),labels=c(1,1.2,1.4,1.6,1.8,2,2.2),
     las=1,cex.axis=1.5,cex.lab=1.4)
mtext("(b) Education Levels",side=3,outer=T,line=-2,cex=1.7)
dev.off()

####

pdf(file="fig1c.pdf")

d=c("High", "Med.", "Low")  
b=c("Unemployment Rate","Mean subjective Unemployment Risk")  
corrs=matrix(nrow=2,ncol=3)
dimnames(corrs)=list(b,d)

corrs[1,1:3]=c(5.49,3.74,3.08)
corrs[2,1]=(mean(a$subj_urisk[which(a$ticino==1|a$geneva==1)]))-1
corrs[2,2]=(mean(a$subj_urisk[which(a$mittell==1|a$nw==1|a$zurich==1)]))-1                  
corrs[2,3]=(mean(a$subj_urisk[which(a$central==1|a$east==1)]))-1

par(mfrow=c(1,2))
barplot(corrs[1,],col=c("gray60","gray60","gray60"),las=1,ylab="Mean Unemployment Rate (Per Cent)",
        ylim=c(0,6),xlab="",
        ,cex.lab=1.4,cex.axis=1.5,cex.names=1.1)
barplot(corrs[2,],col=c("gray60","gray60","gray60"),
        ylab="Mean Subjective Unemployment Risk",ylim=c(0,1.2),yaxt="n"
        ,xlab="",cex.lab=1.4,cex.names=1.1)
axis(2,at=c(0,0.2,0.4,0.6,0.8,1,1.2),labels=c(1,1.2,1.4,1.6,1.8,2,2.2),las=1,cex.axis=1.5)
#mtext("Region (Classified by Unemployment Rate)",side=1,outer=T,line=-3,cex=1.4)
#mtext("(Classified by Unemployment Rate)",side=1,outer=T,line=-2,cex=1.4)
mtext("(c) Regions (Classified by Unemployment Rate)",side=3,outer=T,line=-2,cex=1.7)
dev.off()

## 4. Occupation ##

a=a[!is.na(a$is1maj),]
a$is1maj=as.numeric(a$is1maj)

pdf(file="fig1d.pdf")
d=c("Managers", "Professionals", "Technicians","Clerical","Services & sales"
    ,"Craft & Trades","Manual","Elementary")  
b=c("Unemployment Rate","Mean subjective Unemployment Risk")  
corrs=matrix(nrow=2,ncol=8)
dimnames(corrs)=list(b,d)

corrs[1,1:8]=c(2.86,  1.89,	1.86,	3.54,	4.13,	3.02,	3.66,	3.65)
corrs[2,1]=(mean(a$subj_urisk[which(a$is1maj==6)]))-1
corrs[2,2]=(mean(a$subj_urisk[which(a$is1maj==7)]))-1
corrs[2,3]=(mean(a$subj_urisk[which(a$is1maj==8)]))-1                  
corrs[2,4]=(mean(a$subj_urisk[which(a$is1maj==9)]))-1  
corrs[2,5]=(mean(a$subj_urisk[which(a$is1maj==10)]))-1  
corrs[2,6]=(mean(a$subj_urisk[which(a$is1maj==12)]))-1  
corrs[2,7]=(mean(a$subj_urisk[which(a$is1maj==13)]))-1  
corrs[2,8]=(mean(a$subj_urisk[which(a$is1maj==14)]))-1  

par(mar=c(4,9,3,2))
par(mfrow=c(2,1))
barplot(corrs[1,],horiz=T,col=c("gray60","gray60","gray60",
                                "gray60","gray60","gray60","gray60","gray60"),
        las=1,xlab="Mean Unemployment Rate",
        xlim=c(0,5),ylab="",cex.lab=1.4,cex.axis=1.5,cex.names=1.2)
barplot(corrs[2,],horiz=T,col=c("gray60","gray60","gray60",
                                "gray60","gray60","gray60","gray60","gray60"),las=1,
        xlab="Mean Subjective Unemployment Risk (SHP)",xlim=c(0,1.4),cex.lab=1.4,cex.names=1.2,
        xaxt="n",ylab="")
axis(1,at=c(0,0.2,0.4,0.6,0.8,1,1.2,1.4),labels=c(1,1.2,1.4,1.6,1.8,2,2.2,2.4),cex.axis=1.5)
mtext("(d) Occupations",side=3,outer=T,line=-2,cex=1.7)

#mtext("Sources: Swiss Household Panel (SHP), Eurostat and Swiss Federal Statistical Office",side=1,outer=T,line=-1,cex=0.65)
dev.off()

#####################################################################
#####################################################################



# Figure 2: Unemployment Risk Over Time

pdf(file="fig2.pdf",width=6.5,height=4.5)
par(mar=c(5,4,2,5))
#par(oma = c(4, 1, 3, 1))

x=c(1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2011)
surate=c(mean(a$subj_urisk[a$wave==1]),mean(a$subj_urisk[a$wave==2]),mean(a$subj_urisk[a$wave==3]),
         mean(a$subj_urisk[a$wave==4]),mean(a$subj_urisk[a$wave==5]),mean(a$subj_urisk[a$wave==6]),
         mean(a$subj_urisk[a$wave==7]),mean(a$subj_urisk[a$wave==8]),mean(a$subj_urisk[a$wave==9]),
         mean(a$subj_urisk[a$wave==10]),mean(a$subj_urisk[a$wave==11]),mean(a$subj_urisk[a$wave==12]))
urate=c(3.1,2.7,2.5,2.9,4.1,4.3,4.4,4.0,3.6,3.4,4.1,3.6)

plot(x,surate,type="l",lty=1,lwd=2,ylab="", xlab="",cex.axis=0.8,ylim=c(1.4,2.2))
mtext("Mean Subjective Unemployment Risk",side=2,line=2.5,cex=1)
#mtext("Unemployment Risk",side=2,line=2,cex=0.6)
par(new=T)
plot(x,urate,type="l",lty=2,lwd=2,xaxt="n",yaxt="n",ylab="",xlab="",ylim=c(2.4,4.5))
axis(4)
mtext("Swiss Unemployment Rate (%)",side=4,line=2.5,cex=1)

par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend(-0.9,-0.85, c("National Unemployment Rate (%)", "Mean Subjective Unemployment Risk"), 
       xpd = TRUE, horiz = TRUE,inset = c(0,0), bty = "n", lty=c(1,2),cex = 0.8)
#mtext("Source: Swiss Household Panel (SHP) and Swiss Federal Statistics Office",side=1,outer=T,line=-1,cex=0.7)

dev.off()

#####################################################################
#####################################################################



# Figure 3: Evolution over Time

#(a) aggregate prefs, those who began at wave 1
yrs=c(1999:2009,2011)

ss2=ss0=t2=t0=ss2b=ss0b=t2b=t0b=c()

for(i in 1:12){
  ss2[i]=length(a$opn_socsec[which(a$opn_socsec==1&a$wave==i&a$len1==1)])/length(a$opn_socsec[which(a$wave==i&a$len1==1)])*100
  ss0[i]=length(a$opn_socsec[which(a$opn_socsec==0&a$wave==i&a$len1==1)])/length(a$opn_socsec[which(a$wave==i&a$len1==1)])*100
  t2[i]=length(a$opn_taxes[which(a$opn_taxes==1&a$wave==i&a$len1==1)])/length(a$opn_taxes[which(a$wave==i&a$len1==1)])*100
  t0[i]=length(a$opn_taxes[which(a$opn_taxes==0&a$wave==i&a$len1==1)])/length(a$opn_taxes[which(a$wave==i&a$len1==1)])*100
  ss2b[i]=length(a$opn_socsec[which(a$opn_socsec==1&a$wave==i&a$len1==1&a$nobs==12)])/length(a$opn_socsec[which(a$wave==i&a$len1==1&a$nobs==12)])*100
  ss0b[i]=length(a$opn_socsec[which(a$opn_socsec==0&a$wave==i&a$len1==1&a$nobs==12)])/length(a$opn_socsec[which(a$wave==i&a$len1==1&a$nobs==12)])*100
  t2b[i]=length(a$opn_taxes[which(a$opn_taxes==1&a$wave==i&a$len1==1&a$nobs==12)])/length(a$opn_taxes[which(a$wave==i&a$len1==1&a$nobs==12)])*100
  t0b[i]=length(a$opn_taxes[which(a$opn_taxes==0&a$wave==i&a$len1==1&a$nobs==12)])/length(a$opn_taxes[which(a$wave==i&a$len1==1&a$nobs==12)])*100
}

pdf(file="fig3.pdf")
par(mfrow=c(2,2),mar=c(7.5,4,4,2))

plot(yrs,ss2,type="o",lwd=1.5,pch=16,col="gray25",ylim=c(-5,100), xlab="",
     cex.axis=0.8,yaxt="n",
     main="(a) Aggregate Preferences over Time, \n those who started at Wave 1",
     cex.main=0.95,font.main=1,ylab="Percentage")
lines(yrs,ss2b,type="o",lty=2,pch=16,lwd=1.5,col="gray25")
lines(yrs,ss0,type="o",pch=15,col="gray25")
lines(yrs,ss0b,type="o",lty=2,pch=15,lwd=1.5,col="gray25")
lines(yrs,t2,type="o",pch=17,col="gray25")
lines(yrs,t2b,type="o",lty=2,pch=17,lwd=1.5,col="gray25")
lines(yrs,t0,type="o",pch=18,col="gray25")
lines(yrs,t0b,type="o",lty=2,pch=18,lwd=1.5,col="gray25")
axis(2, las=1)


#############

#(b) aggregate prefs, those who began at wave 6
yrs=c(2004:2009,2011)

ss2=ss0=t2=t0=ss2b=ss0b=t2b=t0b=c()

for(i in 1:7){
  ss2[i]=length(a$opn_socsec[which(a$opn_socsec==1&a$wave==(i+5)&a$len1==6)])/length(a$opn_socsec[which(a$wave==(i+5)&a$len1==6)])*100
  ss0[i]=length(a$opn_socsec[which(a$opn_socsec==0&a$wave==(i+5)&a$len1==6)])/length(a$opn_socsec[which(a$wave==(i+5)&a$len1==6)])*100
  t2[i]=length(a$opn_taxes[which(a$opn_taxes==1&a$wave==(i+5)&a$len1==6)])/length(a$opn_taxes[which(a$wave==(i+5)&a$len1==6)])*100
  t0[i]=length(a$opn_taxes[which(a$opn_taxes==0&a$wave==(i+5)&a$len1==6)])/length(a$opn_taxes[which(a$wave==(i+5)&a$len1==6)])*100
  ss2b[i]=length(a$opn_socsec[which(a$opn_socsec==1&a$wave==(i+5)&a$len1==6&a$nobs==7)])/length(a$opn_socsec[which(a$wave==(i+5)&a$len1==6&a$nobs==7)])*100
  ss0b[i]=length(a$opn_socsec[which(a$opn_socsec==0&a$wave==(i+5)&a$len1==6&a$nobs==7)])/length(a$opn_socsec[which(a$wave==(i+5)&a$len1==6&a$nobs==7)])*100
  t2b[i]=length(a$opn_taxes[which(a$opn_taxes==1&a$wave==(i+5)&a$len1==6&a$nobs==7)])/length(a$opn_taxes[which(a$wave==(i+5)&a$len1==6&a$nobs==7)])*100
  t0b[i]=length(a$opn_taxes[which(a$opn_taxes==0&a$wave==(i+5)&a$len1==6&a$nobs==7)])/length(a$opn_taxes[which(a$wave==(i+5)&a$len1==6&a$nobs==7)])*100
}

plot(yrs,ss2,type="o",lwd=1.5,pch=16,col="gray25",ylim=c(-5,100), xlab="",
     cex.axis=0.75,yaxt="n",
     main="(b) Aggregate Preferences over Time, \n those who started at Wave 6"
     ,cex.main=0.95,font.main=1,ylab="Percentage")
lines(yrs,ss2b,type="o",lty=2,pch=16,lwd=1.5,col="gray25")
lines(yrs,ss0,type="o",pch=15,col="gray25")
lines(yrs,ss0b,type="o",lty=2,pch=15,lwd=1.5,col="gray25")
lines(yrs,t2,type="o",pch=17,col="gray25")
lines(yrs,t2b,type="o",lty=2,pch=17,lwd=1.5,col="gray25")
lines(yrs,t0,type="o",pch=18,col="gray25")
lines(yrs,t0b,type="o",lty=2,pch=18,lwd=1.5,col="gray25")
axis(2, las=1)

###########

# Bottom two panels - stability

# 10 year period

yrs=c(1,2,3,4,5,6,7,8,9,10,12)

ch0=ch0t=ch0b=ch0bt=ch1=ch1t=ch1b=ch1bt=c()

for (i in 1:11){
  ch0[i]=(length(a$opn_socsec[which(a$wavegap==i&a$beliefdiff==0&a$len1==1)])/length(a$opn_socsec[which(a$wavegap==i&a$len1==1)]))*100
  ch0t[i]=(length(a$opn_taxes[which(a$wavegap==i&a$beliefdifft==0&a$len1==1)])/length(a$opn_taxes[which(a$wavegap==i&a$len1==1)]))*100
  ch0b[i]=(length(a$opn_socsec[which(a$wavegap==i&a$beliefdiff==0&a$len1==1&a$nobs==12)])/length(a$opn_socsec[which(a$wavegap==i&a$len1==1&a$nobs==12)]))*100
  ch0bt[i]=(length(a$opn_taxes[which(a$wavegap==i&a$beliefdifft==0&a$len1==1&a$nobs==12)])/length(a$opn_taxes[which(a$wavegap==i&a$len1==1&a$nobs==12)]))*100
  ch1[i]=(length(a$opn_socsec[which(a$wavegap==i&(a$beliefdiff==0.5|a$beliefdiff==-0.5)&a$len1==1)])/length(a$opn_socsec[which(a$wavegap==i&a$len1==1)]))*100
  ch1t[i]=(length(a$opn_taxes[which(a$wavegap==i&(a$beliefdifft==0.5|a$beliefdifft==-0.5)&a$len1==1)])/length(a$opn_taxes[which(a$wavegap==i&a$len1==1)]))*100
  ch1b[i]=(length(a$opn_socsec[which(a$wavegap==i&(a$beliefdiff==0.5|a$beliefdiff==-0.5)&a$len1==1&a$nobs==12)])/length(a$opn_socsec[which(a$wavegap==i&a$len1==1&a$nobs==12)]))*100
  ch1bt[i]=(length(a$opn_taxes[which(a$wavegap==i&(a$beliefdifft==0.5|a$beliefdifft==-0.5)&a$len1==1&a$nobs==12)])/length(a$opn_taxes[which(a$wavegap==i&a$len1==1&a$nobs==12)]))*100
}

plot(yrs,ch0,type="o",lwd=1.5,yaxt="n",pch=16,col="gray25",ylim=c(0,100), 
     xlab="",ylab="Percent",
     main="(c) Stability of Preferences over Time, \n those who started at Wave 1",
     cex.main=0.95,font.main=1,cex.lab=0.85)
lines(yrs,ch0b,type="o",lty=2,pch=16,lwd=1.5,col="gray25")
lines(yrs,ch1,type="o",pch=15,lwd=1.5,col="gray25")
lines(yrs,ch1b,type="o",pch=15,lwd=1.5,lty=2,col="gray25")
lines(yrs,ch0t,type="o",pch=17,lwd=1.5,col="gray25")
lines(yrs,ch0bt,type="o",pch=17,lwd=1.5,lty=2,col="gray25")
lines(yrs,ch1t,type="o",pch=18,lwd=1.5,col="gray25")
lines(yrs,ch1bt,type="o",pch=18,lwd=1.5,lty=2,col="gray25")
axis(2, las=1)
title(xlab="Years since first appearance in the panel", line=2, cex.lab=0.9)

#############

# 6 year period

yrs=c(1,2,3,4,5,7)
ch0=ch0t=ch0b=ch0bt=ch1=ch1t=ch1b=ch1bt=c()

for (i in 1:6){
  ch0[i]=(length(a$opn_socsec[which(a$wavegap==i&a$beliefdiff==0&a$len1==6)])/length(a$opn_socsec[which(a$wavegap==i&a$len1==6)]))*100
  ch0t[i]=(length(a$opn_taxes[which(a$wavegap==i&a$beliefdifft==0&a$len1==6)])/length(a$opn_taxes[which(a$wavegap==i&a$len1==6)]))*100
  ch0b[i]=(length(a$opn_socsec[which(a$wavegap==i&a$beliefdiff==0&a$len1==6&a$nobs==7)])/length(a$opn_socsec[which(a$wavegap==i&a$len1==6&a$nobs==7)]))*100
  ch0bt[i]=(length(a$opn_taxes[which(a$wavegap==i&a$beliefdifft==0&a$len1==6&a$nobs==7)])/length(a$opn_taxes[which(a$wavegap==i&a$len1==6&a$nobs==7)]))*100
  ch1[i]=(length(a$opn_socsec[which(a$wavegap==i&(a$beliefdiff==0.5|a$beliefdiff==-0.5)&a$len1==6)])/length(a$opn_socsec[which(a$wavegap==i&a$len1==6)]))*100
  ch1t[i]=(length(a$opn_taxes[which(a$wavegap==i&(a$beliefdifft==0.5|a$beliefdifft==-0.5)&a$len1==6)])/length(a$opn_taxes[which(a$wavegap==i&a$len1==6)]))*100
  ch1b[i]=(length(a$opn_socsec[which(a$wavegap==i&(a$beliefdiff==0.5|a$beliefdiff==-0.5)&a$len1==6&a$nobs==7)])/length(a$opn_socsec[which(a$wavegap==i&a$len1==6&a$nobs==7)]))*100
  ch1bt[i]=(length(a$opn_taxes[which(a$wavegap==i&(a$beliefdifft==0.5|a$beliefdifft==-0.5)&a$len1==6&a$nobs==7)])/length(a$opn_taxes[which(a$wavegap==i&a$len1==6&a$nobs==7)]))*100
}


plot(yrs,ch0,type="o",lwd=1.5,yaxt="n",pch=16,col="gray25",ylim=c(0,100), 
     xlab="",
     ylab="Percent",main="(d) Stability of Preferences over Time, \n those who started at Wave 6",
     cex.main=0.95,font.main=1,cex.lab=0.85)
lines(yrs,ch0b,type="o",lty=2,pch=16,lwd=1.5,col="gray25")
lines(yrs,ch1,type="o",pch=15,lwd=1.5,col="gray25")
lines(yrs,ch1b,type="o",pch=15,lwd=1.5,lty=2,col="gray25")
lines(yrs,ch0t,type="o",pch=17,lwd=1.5,col="gray25")
lines(yrs,ch0bt,type="o",pch=17,lwd=1.5,lty=2,col="gray25")
lines(yrs,ch1t,type="o",pch=18,lwd=1.5,col="gray25")
lines(yrs,ch1bt,type="o",pch=18,lwd=1.5,lty=2,col="gray25")
axis(2, las=1)
title(xlab="Years since first appearance in the panel", line=2, cex.lab=0.9)


par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend(-0.5,0.23, c("Higher taxes on rich: support", "Higher social spending: support",
                    "Higher social spending: oppose","Higher taxes on rich: oppose"), 
       xpd = TRUE, inset = c(0,0), lty=c(1,1,1,1),
       pch=c(17,15,16,18),cex = 0.85,ncol=2)

par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend(-0.5,-0.92, c("Taxes on rich: no change", "Social spending: 1 pt. change",
                     "Social spending: no change","Taxes on rich: 1 pt. change"), 
       xpd = TRUE, inset = c(0,0), lty=c(1,1,1,1),
       pch=c(17,15,16,18),cex = 0.85,ncol=2)
dev.off()



#####################################################################
#####################################################################



# Figure 4: Big Changes Fig

a$inc_chdcat[a$inc_chdcat==0]<-NA
a$inc_chupcat[a$inc_chupcat==0]<-NA

##########

# 1. Soc sec and rises in urisk 

#0-1
baseline <- (length(a$diff_socsec[which(a$diff_urisk<2&a$diff_urisk>(-1)&a$socsec_lag!=1&a$diff_socsec>0)])/
               length(a$diff_socsec[which(a$diff_urisk<2&a$diff_urisk>(-1)&a$socsec_lag!=1)]))*100

b1=( (length(a$diff_socsec[which(a$diff_urisk>5&a$socsec_lag!=1&a$diff_socsec>0)])/
        length(a$diff_socsec[which(a$diff_urisk>5&a$socsec_lag!=1)]))*100 ) - baseline   #5+
b2=( (length(a$diff_socsec[which(a$diff_urisk<6&a$diff_urisk>3&a$socsec_lag!=1&a$diff_socsec>0)])/
        length(a$diff_socsec[which(a$diff_urisk<6&a$diff_urisk>3&a$socsec_lag!=1)]))*100 ) - baseline #4-5
b3=( (length(a$diff_socsec[which(a$diff_urisk<4&a$diff_urisk>1&a$socsec_lag!=1&a$diff_socsec>0)])/  #2-3
        length(a$diff_socsec[which(a$diff_urisk<4&a$diff_urisk>1&a$socsec_lag!=1)]))*100 ) - baseline


pdf(file="fig4.pdf",width=8,height=5.5)
par(mfrow=c(2,4))
par(oma = c(4, 0, 3, 1),mar=c(4,4,4,1))

change=c(b3,b2,b1)
plot(change,type="o",ylim=c(-2,8),pch=16,
     main="Rises in Unemployment Risk \n and Social Spending",cex.main=0.9,las=2,cex.axis=0.9,
     xlab="",xaxt="n",ylab="")
abline(h=0,lty=3)
axis(1,at=c(1,2,3),labels=c("2 - 3","4 -5","6 - 10"),cex.axis=0.9)
title(ylab="Percent", line=2, cex.lab=0.9)

###

# 2. Soc sec and falls in urisk 

#0- (-1)
baseline <- (length(a$diff_socsec[which(a$diff_urisk>(-2)&a$diff_urisk<1&a$socsec_lag!=0&a$diff_socsec<0)])/
               length(a$diff_socsec[which(a$diff_urisk>(-2)&a$diff_urisk<1&a$socsec_lag!=0)]))*100

d1=( (length(a$diff_socsec[which(a$diff_urisk<(-5)&a$socsec_lag!=0&a$diff_socsec<0)])/
        length(a$diff_socsec[which(a$diff_urisk<(-5)&a$socsec_lag!=0)]))*100 ) - baseline  #6+
d2=( (length(a$diff_socsec[which(a$diff_urisk>(-6)&a$diff_urisk<(-3)&a$socsec_lag!=0&a$diff_socsec<0)])/
        length(a$diff_socsec[which(a$diff_urisk>(-6)&a$diff_urisk<(-3)&a$socsec_lag!=0)]))*100 ) - baseline #4-5
d3=( (length(a$diff_socsec[which(a$diff_urisk>(-4)&a$diff_urisk<(-1)&a$socsec_lag!=0&a$diff_socsec<0)])/ #2-3
        length(a$diff_socsec[which(a$diff_urisk>(-4)&a$diff_urisk<(-1)&a$socsec_lag!=0)]))*100 ) - baseline 

change=c(d1,d2,d3)
plot(change,type="o",ylim=c(-2,8),pch=17,lty=2,
     main="Falls in Unemployment Risk \n and Social Spending",cex.main=0.9,las=2,cex.axis=0.9,
     xlab="",xaxt="n",ylab="")
abline(h=0,lty=3)
axis(1,at=c(1,2,3),labels=c("-6 to -10","-5 to -4","-3 to -2"),cex.axis=0.9)


############

# 3. Taxes and rises in urisk 

baseline <- (length(a$diff_taxes[which(a$diff_urisk<2&a$diff_urisk>(-1)&a$taxes_lag!=1&a$diff_taxes>0)])/
               length(a$diff_taxes[which(a$diff_urisk<2&a$diff_urisk>(-1)&a$taxes_lag!=1)]))*100

e1=( (length(a$diff_taxes[which(a$diff_urisk>5&a$taxes_lag!=1&a$diff_taxes>0)])/
        length(a$diff_taxes[which(a$diff_urisk>5&a$taxes_lag!=1)]))*100 ) - baseline  #6+
e2=( (length(a$diff_taxes[which(a$diff_urisk<6&a$diff_urisk>3&a$taxes_lag!=1&a$diff_taxes>0)])/
        length(a$diff_taxes[which(a$diff_urisk<6&a$diff_urisk>3&a$taxes_lag!=1)]))*100 ) - baseline #4-5
e3=( (length(a$diff_taxes[which(a$diff_urisk<4&a$diff_urisk>1&a$taxes_lag!=1&a$diff_taxes>0)])/ #2-3
        length(a$diff_taxes[which(a$diff_urisk<4&a$diff_urisk>1&a$taxes_lag!=1)]))*100 ) - baseline


change=c(e3,e2,e1)
plot(change,type="o",ylim=c(-2,8),pch=16,
     main="Rises in Unemployment Risk \n and Taxing the Rich",cex.main=0.9,las=2,cex.axis=0.9,
     xlab="",xaxt="n",ylab="")
abline(h=0,lty=3)
axis(1,at=c(1,2,3),labels=c("2 - 3","4 -5","6 - 10"),cex.axis=0.9)


###

# 4. Taxes and falls in urisk

baseline <- (length(a$diff_taxes[which(a$diff_urisk>(-2)&a$diff_urisk<1&a$taxes_lag!=0&a$diff_taxes<0)])/
               length(a$diff_taxes[which(a$diff_urisk>(-2)&a$diff_urisk<1&a$taxes_lag!=0)]))*100


f1=( (length(a$diff_taxes[which(a$diff_urisk<(-5)&a$taxes_lag!=0&a$diff_taxes<0)])/
        length(a$diff_taxes[which(a$diff_urisk<(-5)&a$taxes_lag!=0)]))*100 ) - baseline  #6+
f2=( (length(a$diff_taxes[which(a$diff_urisk>(-6)&a$diff_urisk<(-3)&a$taxes_lag!=0&a$diff_taxes<0)])/
        length(a$diff_taxes[which(a$diff_urisk>(-6)&a$diff_urisk<(-3)&a$taxes_lag!=0)]))*100 ) - baseline #4-5
f3=( (length(a$diff_taxes[which(a$diff_urisk>(-4)&a$diff_urisk<(-1)&a$taxes_lag!=0&a$diff_taxes<0)])/ #2-3
        length(a$diff_taxes[which(a$diff_urisk>(-4)&a$diff_urisk<(-1)&a$taxes_lag!=0)]))*100 ) - baseline

change=c(f1,f2,f3)
plot(change,type="o",ylim=c(-2,8),pch=17,lty=2,
     main="Falls in Unemployment Risk \n and Taxing the Rich",cex.main=0.9,las=2,cex.axis=0.9,
     xlab="",xaxt="n",ylab="")
abline(h=0,lty=3)
axis(1,at=c(1,2,3),labels=c("-6 to -10","-5 to -4","-3 to -2"),cex.axis=0.9)


#########################


# 5. Soc sec and falls in income 

baseline = (length(a$diff_socsec[which(a$inc_chdcat>(7)&a$socsec_lag!=1&a$diff_socsec>0)])/
              length(a$diff_socsec[which(a$inc_chdcat>(7)&a$socsec_lag!=1)]))*100 #1-20


g1=( (length(a$diff_socsec[which(a$inc_chdcat==1&a$socsec_lag!=1&a$diff_socsec>0)])/
        length(a$diff_socsec[which(a$inc_chdcat==1&a$socsec_lag!=1)]))*100) - baseline #lt 50
g2=( (length(a$diff_socsec[which(a$inc_chdcat>(1)&a$inc_chdcat<(4)&a$socsec_lag!=1&a$diff_socsec>0)])/
        length(a$diff_socsec[which(a$inc_chdcat>(1)&a$inc_chdcat<(4)&a$socsec_lag!=1)]))*100) - baseline #40-50
g3=( (length(a$diff_socsec[which(a$inc_chdcat>(3)&a$inc_chdcat<(6)&a$socsec_lag!=1&a$diff_socsec>0)])/
        length(a$diff_socsec[which(a$inc_chdcat>(3)&a$inc_chdcat<(6)&a$socsec_lag!=1)]))*100) - baseline #30-40
g4=( (length(a$diff_socsec[which(a$inc_chdcat>(5)&a$inc_chdcat<(8)&a$socsec_lag!=1&a$diff_socsec>0)])/
        length(a$diff_socsec[which(a$inc_chdcat>(5)&a$inc_chdcat<(8)&a$socsec_lag!=1)]))*100) - baseline #20-30


change=c(g1,g2,g3,g4)
plot(change,type="o",ylim=c(-2,8),pch=16,
     main="Falls in Income \n and Social Spending",cex.main=0.9,las=2,cex.axis=0.9,
     xlab="",xaxt="n",ylab="")
axis(1,at=c(1,2,3,4),labels=F)
text(c(1,2,3,4), par("usr")[3] - 2, labels = c("<(-50)%","-50% to -40%","-40% to -30%","-30% to -20%")
     , srt = 50, pos = 1, xpd = TRUE,cex=0.8)
abline(h=0,lty=3)
title(ylab="Percent", line=2, cex.lab=0.9)

######

# 6. socsec and rises in income 

baseline <- (length(a$diff_socsec[which(a$inc_chupcat<(5)&a$socsec_lag!=0&a$diff_socsec<0)])/
               length(a$diff_socsec[which(a$inc_chupcat<(5)&a$socsec_lag!=0)]))*100 #<20

h1=( (length(a$diff_socsec[which(a$inc_chupcat==11&a$socsec_lag!=0&a$diff_socsec<0)])/
        length(a$diff_socsec[which(a$inc_chupcat==11&a$socsec_lag!=0)]))*100) - baseline #50+
h2=( (length(a$diff_socsec[which(a$inc_chupcat>(8)&a$inc_chupcat<(11)&a$socsec_lag!=0&a$diff_socsec<0)])/
        length(a$diff_socsec[which(a$inc_chupcat>(8)&a$inc_chupcat<(11)&a$socsec_lag!=0)]))*100) - baseline  #40-50
h3=( (length(a$diff_socsec[which(a$inc_chupcat>(6)&a$inc_chupcat<(9)&a$socsec_lag!=0&a$diff_socsec<0)])/
        length(a$diff_socsec[which(a$inc_chupcat>(6)&a$inc_chupcat<(9)&a$socsec_lag!=0)]))*100) - baseline  #30-40
h4=( (length(a$diff_socsec[which(a$inc_chupcat>(4)&a$inc_chupcat<(7)&a$socsec_lag!=0&a$diff_socsec<0)])/
        length(a$diff_socsec[which(a$inc_chupcat>(4)&a$inc_chupcat<(7)&a$socsec_lag!=0)]))*100) - baseline #20-30


change=c(h4,h3,h2,h1)
plot(change,type="o",ylim=c(-2,8),pch=17,lty=2,
     main="Rises in Income \n and Social Spending",cex.main=0.9,las=2,cex.axis=0.9,
     xlab="",xaxt="n",ylab="")
axis(1,at=c(1,2,3,4),labels=F)
text(c(1,2,3,4), par("usr")[3] - 1.8, labels = c("20% to 30%","30% to 40%","40% to 50%",">50%")
     , srt = 50, pos = 1, xpd = TRUE,cex=0.8)
abline(h=0,lty=3)

###

# 7. taxes and falls in income 

baseline <- (length(a$diff_taxes[which(a$inc_chdcat>(7)&a$taxes_lag!=1&a$diff_taxes>0)])/
               length(a$diff_taxes[which(a$inc_chdcat>(7)&a$taxes_lag!=1)]))*100 #1-20

j1=( (length(a$diff_taxes[which(a$inc_chdcat==1&a$taxes_lag!=1&a$diff_taxes>0)])/
        length(a$diff_taxes[which(a$inc_chdcat==1&a$taxes_lag!=1)]))*100 ) - baseline #lt 50
j2=( (length(a$diff_taxes[which(a$inc_chdcat>(1)&a$inc_chdcat<(4)&a$taxes_lag!=1&a$diff_taxes>0)])/
        length(a$diff_taxes[which(a$inc_chdcat>(1)&a$inc_chdcat<(4)&a$taxes_lag!=1)]))*100 ) - baseline #40-50
j3=( (length(a$diff_taxes[which(a$inc_chdcat>(3)&a$inc_chdcat<(6)&a$taxes_lag!=1&a$diff_taxes>0)])/
        length(a$diff_taxes[which(a$inc_chdcat>(3)&a$inc_chdcat<(6)&a$taxes_lag!=1)]))*100) - baseline #30-40
j4=( (length(a$diff_taxes[which(a$inc_chdcat>(5)&a$inc_chdcat<(8)&a$taxes_lag!=1&a$diff_taxes>0)])/
        length(a$diff_taxes[which(a$inc_chdcat>(5)&a$inc_chdcat<(8)&a$taxes_lag!=1)]))*100 ) - baseline #20-30

change=c(j1,j2,j3,j4)
plot(change,type="o",ylim=c(-2,10),pch=16,
     main="Falls in Income \n and Taxing the Rich",cex.main=0.9,las=2,cex.axis=0.9,
     xlab="",xaxt="n",ylab="")
axis(1,at=c(1,2,3,4),labels=F)
text(c(1,2,3,4), par("usr")[3] - 2, labels = c("<(-50)%","-50% to -40%","-40% to -30%","-30% to -20%")
     , srt = 50, pos = 1, xpd = TRUE,cex=0.8)
abline(h=0,lty=3)

###

# 8. taxes and rises in income

baseline=(length(a$diff_taxes[which(a$inc_chupcat<(5)&a$taxes_lag!=0&a$diff_taxes<0)])/
            length(a$diff_taxes[which(a$inc_chupcat<(5)&a$taxes_lag!=0)]))*100

k1=( (length(a$diff_taxes[which(a$inc_chupcat==11&a$taxes_lag!=0&a$diff_taxes<0)])/
        length(a$diff_taxes[which(a$inc_chupcat==11&a$taxes_lag!=0)]))*100) - baseline #>50
k2=( (length(a$diff_taxes[which(a$inc_chupcat>(8)&a$inc_chupcat<(11)&a$taxes_lag!=0&a$diff_taxes<0)])/
        length(a$diff_taxes[which(a$inc_chupcat>(8)&a$inc_chupcat<(11)&a$taxes_lag!=0)]))*100) - baseline
k3=( (length(a$diff_taxes[which(a$inc_chupcat>(6)&a$inc_chupcat<(9)&a$taxes_lag!=0&a$diff_taxes<0)])/
        length(a$diff_taxes[which(a$inc_chupcat>(6)&a$inc_chupcat<(9)&a$taxes_lag!=0)]))*100) - baseline
k4=( (length(a$diff_taxes[which(a$inc_chupcat>(4)&a$inc_chupcat<(7)&a$taxes_lag!=0&a$diff_taxes<0)])/
        length(a$diff_taxes[which(a$inc_chupcat>(4)&a$inc_chupcat<(7)&a$taxes_lag!=0)]))*100) - baseline


change=c(k4,k3,k2,k1)
plot(change,type="o",ylim=c(-2,8),pch=17,lty=2,
     main="Rises in Income \n and Taxing the Rich",cex.main=0.9,las=2,cex.axis=0.9,
     xlab="",xaxt="n",ylab="")
axis(1,at=c(1,2,3),labels=F)
text(c(1,2,3,4), par("usr")[3] - 1.8, labels = c("20% to 30%","30% to 40%","40% to 50%",">50%")
     , srt = 50, pos = 1, xpd = TRUE,cex=0.8)
abline(h=0,lty=3)

par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend(-0.3,-0.9, c("Higher Preference", "Lower Preference"), xpd = TRUE, horiz = F,
       inset = c(0,0), pch = c(16,17), lty=c(1,2), cex = 1, ncol=2)
#mtext("Source: Swiss Household Panel (SHP)",side=1,outer=T,line=-1,cex=0.6,adj=0)
dev.off()


#####################################################################
#####################################################################



# Figure 5: Big Changes Regression [results from stata file]

## Soc Spending ##
out1 <- data.frame(matrix(NA,4,3))
colnames(out1) <- c("estimate", "ci_lo", "ci_hi")
out1[4,] <- c(0.00797,-0.00306,0.019)
out1[3,] <- c(-0.0009979,-0.0120506,0.0100549)
out1[2,] <- c(0.0008845,-0.0100985,0.010518)
out1[1,] <- c(-0.001161,-0.01284,0.010518)
out1$num <- c(1:4)
out1$names <- c("Rise in\n Income","Fall in Risk of\n Unemployment",
                "Fall in\n Income","Rise in Risk of\n Unemployment")

p1 <- (ggplot(out1, aes(y=num, x=estimate)) + 
         geom_errorbarh(aes(xmin=ci_lo, xmax=ci_hi), height=.1, col="black") +
         geom_point(aes(y=num, x=estimate), size=2, col="black", shape=21, fill="black") +
         geom_vline(aes(xintercept=0)) +
         scale_x_continuous(name="", limits=c(-0.025,0.025)) + 
         scale_y_discrete(name="",limits=out1$names) +
         theme(axis.text.y = element_text(size = 8)) +
         theme(axis.text.x = element_text(size = 8)) +
         ggtitle("(a) Social Spending \n") +
         theme(plot.title = element_text(size = 10)))

## Taxes ##
out2 <- data.frame(matrix(NA,4,3))
colnames(out2) <- c("estimate", "ci_lo", "ci_hi")
out2[4,] <- c(0.007117,-0.00204,0.0163)
out2[3,] <- c(0.007894,-0.001513,0.0173)
out2[2,] <- c(-0.0017,-0.01061,0.007195)
out2[1,] <- c(0.00152,-0.007971,0.01102)
out2$num <- c(1:4)
out2$names <- c("Rise in\n Income","Fall in Risk of\n Unemployment",
                "Fall in\n Income","Rise in Risk of\n Unemployment")

p2 <- ggplot(out2, aes(y=num, x=estimate)) + 
  geom_errorbarh(aes(xmin=ci_lo, xmax=ci_hi), height=.1, col="black") +
  geom_point(aes(y=num, x=estimate), size=2, col="black", shape=21, fill="black") +
  geom_vline(aes(xintercept=0)) +
  scale_x_continuous(name="", limits=c(-0.025,0.025)) + 
  scale_y_discrete(name="",limits=out2$names) +
  theme(axis.text.y = element_text(size = 8)) +
  theme(axis.text.x = element_text(size = 8)) +
  ggtitle("(b) Taxing the Rich \n") +
  theme(plot.title = element_text(size = 10))



multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

pdf(file="fig5.pdf",height=3.5,width=7.5)
multiplot(p1, p2,cols=2)
dev.off()


#####################################################################
#####################################################################



# Figure 5: Big Changes Regression [results from stata file]

## Soc Spending ##
out1 <- data.frame(matrix(NA,4,3))
colnames(out1) <- c("estimate", "ci_lo", "ci_hi")
out1[4,] <- c(0.00797,-0.00306,0.019)
out1[3,] <- c(-0.0009979,-0.0120506,0.0100549)
out1[2,] <- c(0.0008845,-0.0100985,0.010518)
out1[1,] <- c(-0.001161,-0.01284,0.010518)
out1$num <- c(1:4)
out1$names <- c("Rise in\n Income","Fall in Risk of\n Unemployment",
                "Fall in\n Income","Rise in Risk of\n Unemployment")

p1 <- (ggplot(out1, aes(y=num, x=estimate)) + 
         geom_errorbarh(aes(xmin=ci_lo, xmax=ci_hi), height=.1, col="black") +
         geom_point(aes(y=num, x=estimate), size=2, col="black", shape=21, fill="black") +
         geom_vline(aes(xintercept=0)) +
         scale_x_continuous(name="", limits=c(-0.025,0.025)) + 
         scale_y_discrete(name="",limits=out1$names) +
         theme(axis.text.y = element_text(size = 8)) +
         theme(axis.text.x = element_text(size = 8)) +
         ggtitle("(a) Social Spending \n") +
         theme(plot.title = element_text(size = 10)))

## Taxes ##
out2 <- data.frame(matrix(NA,4,3))
colnames(out2) <- c("estimate", "ci_lo", "ci_hi")
out2[4,] <- c(0.007117,-0.00204,0.0163)
out2[3,] <- c(0.007894,-0.001513,0.0173)
out2[2,] <- c(-0.0017,-0.01061,0.007195)
out2[1,] <- c(0.00152,-0.007971,0.01102)
out2$num <- c(1:4)
out2$names <- c("Rise in\n Income","Fall in Risk of\n Unemployment",
                "Fall in\n Income","Rise in Risk of\n Unemployment")

p2 <- ggplot(out2, aes(y=num, x=estimate)) + 
  geom_errorbarh(aes(xmin=ci_lo, xmax=ci_hi), height=.1, col="black") +
  geom_point(aes(y=num, x=estimate), size=2, col="black", shape=21, fill="black") +
  geom_vline(aes(xintercept=0)) +
  scale_x_continuous(name="", limits=c(-0.025,0.025)) + 
  scale_y_discrete(name="",limits=out2$names) +
  theme(axis.text.y = element_text(size = 8)) +
  theme(axis.text.x = element_text(size = 8)) +
  ggtitle("(b) Taxing the Rich \n") +
  theme(plot.title = element_text(size = 10))



multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

pdf(file="bigchangeregs.pdf",height=3.5,width=7.5)
multiplot(p1, p2,cols=2)
graphics.off()


#####################################################################
#####################################################################


# Figure 6: Family Background

## Soc Spending ##
out <- data.frame(matrix(NA,3,3))
colnames(out) <- c("estimate", "ci_lo", "ci_hi")

test1s <- t.test(a$opn_socsec[a$lwingdad==1], a$opn_socsec[a$lwingdad==0])
out[3,1]<- test1s$estimate[1] - test1s$estimate[2]
out[3,c(2,3)]<- test1s$conf.int[1:2]

test1s <- t.test(a$opn_socsec[a$lwingmum==1], a$opn_socsec[a$lwingmum==0])
out[2,1]<- test1s$estimate[1] - test1s$estimate[2]
out[2,c(2,3)]<- test1s$conf.int[1:2]

test1s <- t.test(a$opn_socsec[a$poor==1], a$opn_socsec[a$poor==0])
out[1,1]<- test1s$estimate[1] - test1s$estimate[2]
out[1,c(2,3)]<- test1s$conf.int[1:2]

out$num <- c(1:3)
out$names <- c("Poor Family","Leftwing Mother","Leftwing Father")

p1 <- ggplot(out, aes(y=num, x=estimate)) + 
  geom_errorbarh(aes(xmin=ci_lo, xmax=ci_hi), height=.1, col="black") +
  geom_point(aes(y=num, x=estimate), size=2, col="black", shape=21, fill="black") +
  geom_vline(aes(xintercept=0)) +
  scale_x_continuous(name="", limits=c(-0.02,0.18)) + 
  scale_y_discrete(name="",limits=out$names) +
  theme(axis.text.y = element_text(size = 8)) +
  theme(axis.text.x = element_text(size = 8)) +
  ggtitle("(a) Social Spending \n") +
  theme(plot.title = element_text(size = 10))

###

## Taxes ##
out <- data.frame(matrix(NA,3,3))
colnames(out) <- c("estimate", "ci_lo", "ci_hi")

test1s <- t.test(a$opn_taxes[a$lwingdad==1], a$opn_taxes[a$lwingdad==0])
out[3,1]<- test1s$estimate[1] - test1s$estimate[2]
out[3,c(2,3)]<- test1s$conf.int[1:2]

test1s <- t.test(a$opn_taxes[a$lwingmum==1], a$opn_taxes[a$lwingmum==0])
out[2,1]<- test1s$estimate[1] - test1s$estimate[2]
out[2,c(2,3)]<- test1s$conf.int[1:2]

test1s <- t.test(a$opn_taxes[a$poor==1], a$opn_taxes[a$poor==0])
out[1,1]<- test1s$estimate[1] - test1s$estimate[2]
out[1,c(2,3)]<- test1s$conf.int[1:2]

out$num <- c(1:3)
out$names <- c("Poor Family","Leftwing Mother","Leftwing Father")


p2 <- ggplot(out, aes(y=num, x=estimate)) + 
  geom_errorbarh(aes(xmin=ci_lo, xmax=ci_hi), height=.1, col="black") +
  geom_point(aes(y=num, x=estimate), size=2, col="black", shape=21, fill="black") +
  geom_vline(aes(xintercept=0)) +
  scale_x_continuous(name="", limits=c(-0.02,0.11)) + 
  scale_y_discrete(name="",limits=out$names) +
  theme(axis.text.y = element_text(size = 8)) +
  theme(axis.text.x = element_text(size = 8)) +
  ggtitle("(b) Taxing the Rich \n") +
  theme(plot.title = element_text(size = 10))

### Combo
pdf(file="fig6.pdf",height=2.6,width=7.5)
multiplot(p1, p2,cols=2)
graphics.off()

#####################################################################
#####################################################################


# Figure 7: Family Background regressions [output from stata]


## Soc Spending, lwingdad ##
out1 <- data.frame(matrix(NA,2,3))
colnames(out1) <- c("estimate", "ci_lo", "ci_hi")
out1[2,] <- c(0.0831,0.0616,0.105)
out1[1,] <- c(0.0122,-0.00724,0.0316)
out1$num <- c(1:2)
out1$names <- c("Poor Family","Leftwing Father")

pr1 <- ggplot(out1, aes(y=num, x=estimate)) + 
  geom_errorbarh(aes(xmin=ci_lo, xmax=ci_hi), height=.1, col="black") +
  geom_point(aes(y=num, x=estimate), size=2, col="black", shape=21, fill="black") +
  geom_vline(aes(xintercept=0)) +
  scale_x_continuous(name="", limits=c(-0.02,0.11)) + 
  scale_y_discrete(name="",limits=out1$names) +
  theme(axis.text.y = element_text(size = 8)) +
  theme(axis.text.x = element_text(size = 8)) +
  ggtitle("(a) Social Spending, Leftwing Father") +
  theme(plot.title = element_text(size = 10))


## Soc Spending, lwingmum ##
out2 <- data.frame(matrix(NA,2,3))
colnames(out2) <- c("estimate", "ci_lo", "ci_hi")
out2[2,] <- c(0.0805,0.0581,0.103)
out2[1,] <- c(0.01988,-0.0011,0.0409)
out2$num <- c(1:2)
out2$names <- c("Poor Family","Leftwing Mother")

pr2 <- ggplot(out2, aes(y=num, x=estimate)) + 
  geom_errorbarh(aes(xmin=ci_lo, xmax=ci_hi), height=.1, col="black") +
  geom_point(aes(y=num, x=estimate), size=2, col="black", shape=21, fill="black") +
  geom_vline(aes(xintercept=0)) +
  scale_x_continuous(name="", limits=c(-0.02,0.11)) + 
  scale_y_discrete(name="",limits=out2$names) +
  theme(axis.text.y = element_text(size = 8)) +
  theme(axis.text.x = element_text(size = 8)) +
  ggtitle("(b)  Social Spending, Leftwing Mother") +
  theme(plot.title = element_text(size = 10))


## Taxes, lwingdad ##
out3 <- data.frame(matrix(NA,2,3))
colnames(out3) <- c("estimate", "ci_lo", "ci_hi")
out3[2,] <- c(0.0297,0.0127,0.0466)
out3[1,] <- c(0.0201,0.005,0.0352)
out3$num <- c(1:2)
out3$names <- c("Poor Family","Leftwing Father")

pr3 <- ggplot(out3, aes(y=num, x=estimate)) + 
  geom_errorbarh(aes(xmin=ci_lo, xmax=ci_hi), height=.1, col="black") +
  geom_point(aes(y=num, x=estimate), size=2, col="black", shape=21, fill="black") +
  geom_vline(aes(xintercept=0)) +
  scale_x_continuous(name="", limits=c(-0.02,0.11)) + 
  scale_y_discrete(name="",limits=out3$names) +
  theme(axis.text.y = element_text(size = 8)) +
  theme(axis.text.x = element_text(size = 8)) +
  ggtitle("(c) Taxing the Rich, Leftwing Father") +
  theme(plot.title = element_text(size = 10))


## Taxes, lwingmum ##
out4 <- data.frame(matrix(NA,2,3))
colnames(out4) <- c("estimate", "ci_lo", "ci_hi")
out4[2,] <- c(0.0319,0.0139,0.050)
out4[1,] <- c(0.0274,0.0107,0.044)
out4$num <- c(1:2)
out4$names <- c("Poor Family","Leftwing Mother")

pr4 <- ggplot(out4, aes(y=num, x=estimate)) + 
  geom_errorbarh(aes(xmin=ci_lo, xmax=ci_hi), height=.1, col="black") +
  geom_point(aes(y=num, x=estimate), size=2, col="black", shape=21, fill="black") +
  geom_vline(aes(xintercept=0)) +
  scale_x_continuous(name="", limits=c(-0.02,0.11)) + 
  scale_y_discrete(name="",limits=out4$names) +
  theme(axis.text.y = element_text(size = 8)) +
  theme(axis.text.x = element_text(size = 8)) +
  ggtitle("(d) Taxing the Rich, Leftwing Mother") +
  theme(plot.title = element_text(size = 10))


# Fig
pdf(file="errorbars2.pdf",height=4,width=7)
multiplot(pr1, pr2, pr3, pr4,cols=2)
dev.off()



#####################################################################
#####################################################################


# Figure 8: Variability by Age

ch0old=ch0oldt=ch0y=ch0yt=c()
ch1old=ch1oldt=ch1y=ch1yt=c()
ch2old=ch2oldt=ch2y=ch2yt=c()

for (i in 1:8){
  ch0old[i]=(length(a$opn_socsec[which(a$wavegap==i&a$beliefdiff==0&a$age1>25&a$noyears>5)])/
               length(a$opn_socsec[which(a$wavegap==i&a$age1>25&a$noyears>5)]))*100
  ch0oldt[i]=(length(a$opn_taxes[which(a$wavegap==i&a$beliefdifft==0&a$age1>25&a$noyears>5)])/
                length(a$opn_taxes[which(a$wavegap==i&a$age1>25&a$noyears>5)]))*100
  ch0y[i]=(length(a$opn_socsec[which(a$wavegap==i&a$beliefdiff==0&a$age1<26&a$noyears>5)])/
             length(a$opn_socsec[which(a$wavegap==i&a$age1<26&a$noyears>5)]))*100
  ch0yt[i]=(length(a$opn_taxes[which(a$wavegap==i&a$beliefdifft==0&a$age1<26&a$noyears>5)])/
              length(a$opn_taxes[which(a$wavegap==i&a$age1<26&a$noyears>5)]))*100
  
  ch1old[i]=(length(a$opn_socsec[which(a$wavegap==i&(a$beliefdiff==0.5|a$beliefdiff==-0.5)&a$age1>25&a$noyears>5)])/
               length(a$opn_socsec[which(a$wavegap==i&a$age1>25&a$noyears>5)]))*100
  ch1oldt[i]=(length(a$opn_taxes[which(a$wavegap==i&(a$beliefdifft==0.5|a$beliefdifft==-0.5)&a$age1>25&a$noyears>5)])/
                length(a$opn_taxes[which(a$wavegap==i&a$age1>25&a$noyears>5)]))*100
  ch1y[i]=(length(a$opn_socsec[which(a$wavegap==i&(a$beliefdiff==0.5|a$beliefdiff==-0.5)&a$age1<26&a$noyears>5)])/
             length(a$opn_socsec[which(a$wavegap==i&a$age1<26&a$noyears>5)]))*100
  ch1yt[i]=(length(a$opn_taxes[which(a$wavegap==i&(a$beliefdifft==0.5|a$beliefdifft==-0.5)&a$age1<26&a$noyears>5)])/
              length(a$opn_taxes[which(a$wavegap==i&a$age1<26&a$noyears>5)]))*100
  
  ch2old[i]=(length(a$opn_socsec[which(a$wavegap==i&(a$beliefdiff==1|a$beliefdiff==-1)&a$age1>25&a$noyears>5)])/
               length(a$opn_socsec[which(a$wavegap==i&a$age1>25&a$noyears>5)]))*100
  ch2oldt[i]=(length(a$opn_taxes[which(a$wavegap==i&(a$beliefdifft==1|a$beliefdifft==-1)&a$age1>25&a$noyears>5)])/
                length(a$opn_taxes[which(a$wavegap==i&a$age1>25&a$noyears>5)]))*100
  ch2y[i]=(length(a$opn_socsec[which(a$wavegap==i&(a$beliefdiff==1|a$beliefdiff==-1)&a$age1<26&a$noyears>5)])/
             length(a$opn_socsec[which(a$wavegap==i&a$age1<26&a$noyears>5)]))*100
  ch2yt[i]=(length(a$opn_taxes[which(a$wavegap==i&(a$beliefdifft==1|a$beliefdifft==-1)&a$age1<26&a$noyears>5)])/
              length(a$opn_taxes[which(a$wavegap==i&a$age1<26&a$noyears>5)]))*100
}


yrs=c(1:8)

pdf(file="varbyage.pdf",width=7.5,height=5.5)
par(mfrow=c(1,2),mar=c(6,4,4,2))
plot(yrs,ch0old,type="o",lwd=2,las=1,pch=16,col="gray25",ylim=c(0,100), 
     xlab="",
     ylab="Percent",main="(a) Higher Social Spending",
     cex.main=0.8,cex.lab=0.8,cex.axis=0.8,yaxt="n",font.main=1)
lines(yrs,ch0y,type="o",lty=2,pch=16,lwd=2,col="gray25")
lines(yrs,ch1old,type="o",pch=17,lwd=2,col="gray25")
lines(yrs,ch1y,type="o",pch=17,lwd=2,lty=2,col="gray25")
lines(yrs,ch2old,type="o",pch=15,lwd=2,col="gray25")
lines(yrs,ch2y,type="o",pch=15,lwd=2,lty=2,col="gray25")
axis(2, las=1,cex.axis=0.8)
title(xlab="Years since first appearance in the panel", line=2, cex.lab=0.9)
text(4.7,66, "No Change", cex=0.6)
text(5,37, "1 pt. change", cex=0.6)
text(2.3,14.5, "2 pt. change", cex=0.6)

plot(yrs,ch0oldt,type="o",lwd=2,las=1,pch=16,col="gray25",ylim=c(0,100), 
     xlab="",
     ylab="Percent",main="(b) Higher Taxes on the Rich",
     cex.main=0.8,cex.lab=0.8,cex.axis=0.8,yaxt="n",font.main=1)
lines(yrs,ch0yt,type="o",lty=2,pch=16,lwd=2,col="gray25")
lines(yrs,ch1oldt,type="o",pch=17,lwd=2,col="gray25")
lines(yrs,ch1yt,type="o",pch=17,lwd=2,lty=2,col="gray25")
lines(yrs,ch2oldt,type="o",pch=15,lwd=2,col="gray25")
lines(yrs,ch2yt,type="o",pch=15,lwd=2,lty=2,col="gray25")
axis(2, las=1,cex.axis=0.8)
title(xlab="Years since first appearance in the panel", line=2, cex.lab=0.9)
text(5,80.5, "No change", cex=0.6)
text(2,24.2, "1 pt. change", cex=0.6)
text(4.5,1.6, "2 pt. change", cex=0.6)

par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend(-0.55,-0.95, c("Respondents Aged 26+", "Respondents Aged 18-25"), 
       xpd = TRUE, horiz = TRUE,inset = c(0,0), lty=c(1,2),cex = 0.8)
dev.off()


#####################################################################
#####################################################################


# Figure 9: Variability by Age II


### soc spending, income ###
# means, old
a1=mean(a$chincome_age[which(a$chbelief==1&a$age1>25&a$yearstodate>5)])
a2=mean(a$chincome_age[which(a$chbelief==0.5&a$age1>25&a$yearstodate>5)])
a3=mean(a$chincome_age[which(a$chbelief==0&a$age1>25&a$yearstodate>5)])
a4=mean(a$chincome_age[which(a$chbelief==-0.5&a$age1>25&a$yearstodate>5)])
a5=mean(a$chincome_age[which(a$chbelief==-1&a$age1>25&a$yearstodate>5)])
meano=c(a5,a4,a3,a2,a1)
totalmedo=mean(a$chincome_age[which(a$age1>25&a$yearstodate>5)])

b1=mean(a$chincome_age[which(a$chbelief==1&a$age1<26&a$yearstodate>5)])
b2=mean(a$chincome_age[which(a$chbelief==0.5&a$age1<26&a$yearstodate>5)])
b3=mean(a$chincome_age[which(a$chbelief==0&a$age1<26&a$yearstodate>5)])
b4=mean(a$chincome_age[which(a$chbelief==-0.5&a$age1<26&a$yearstodate>5)])
b5=mean(a$chincome_age[which(a$chbelief==-1&a$age1<26&a$yearstodate>5)])
meany=c(b5,b4,b3,b2,b1)
totalmedy=mean(a$chincome_age[which(a$age1<26&a$yearstodate>5)])


pdf(file="changesbyage.pdf")
par(oma = c(4, 2, 3, 1))
par(mar=c(5,5,4,2))
par(mfrow=c(2,2))

plot(meano,ylim=c(0,150),type="l",lwd=2,
     xlab="Change in Preference for Higher Social \n Spending since first wave",
     ylab="Mean Percentage Change in \n Income since first wave",
     main="(a) Higher Social Spending and \n Income"
     ,cex.axis=0.9,cex.main=0.9,xaxt="n",cex.lab=0.8,las=2,font.main=1)
box()
axis(1,at=c(1,2,3,4,5),labels=c(-2,-1,0,1,2),cex.axis=0.9)
lines(meany,type="l",lty=2,lwd=2)


### taxes, income ###

a1=mean(a$chincome_age[which(a$chbelieft==1&a$age1>25&a$yearstodate>5)])
a2=mean(a$chincome_age[which(a$chbelieft==0.5&a$age1>25&a$yearstodate>5)])
a3=mean(a$chincome_age[which(a$chbelieft==0&a$age1>25&a$yearstodate>5)])
a4=mean(a$chincome_age[which(a$chbelieft==-0.5&a$age1>25&a$yearstodate>5)])
a5=mean(a$chincome_age[which(a$chbelieft==-1&a$age1>25&a$yearstodate>5)])
meano=c(a5,a4,a3,a2,a1)
totalmedo=mean(a$chincome_age[which(a$age1>25&a$yearstodate>5)])

b1=mean(a$chincome_age[which(a$chbelieft==1&a$age1<26&a$yearstodate>5)])
b2=mean(a$chincome_age[which(a$chbelieft==0.5&a$age1<26&a$yearstodate>5)])
b3=mean(a$chincome_age[which(a$chbelieft==0&a$age1<26&a$yearstodate>5)])
b4=mean(a$chincome_age[which(a$chbelieft==-0.5&a$age1<26&a$yearstodate>5)])
b5=mean(a$chincome_age[which(a$chbelieft==-1&a$age1<26&a$yearstodate>5)])
meany=c(b5,b4,b3,b2,b1)
totalmedy=mean(a$chincome_age[which(a$age1<26&a$yearstodate>5)])

plot(meano,ylim=c(0,200),type="l",lwd=2,
     xlab="Change in Preference for Higher Taxes \n on Rich since first wave",
     ylab="Mean Percentage Change in \n Income since first wave",
     main="(b) Higher Taxes on the Rich and \n Income",
     cex.axis=0.9,cex.main=0.9,xaxt="n",cex.lab=0.8,las=2,font.main=1)
box()
axis(1,at=c(1,2,3,4,5),labels=c(-2,-1,0,1,2),cex.axis=0.9)
lines(meany,type="l",lty=2,lwd=2)
#abline(h=0)

######################################################################

#U Risk

### Soc spending, urisk ###

# Mens, old
a1=mean(a$churisk_age[which(a$chbelief==1&a$age1>25&a$yearstodate>5)])
a2=mean(a$churisk_age[which(a$chbelief==0.5&a$age1>25&a$yearstodate>5)])
a3=mean(a$churisk_age[which(a$chbelief==0&a$age1>25&a$yearstodate>5)])
a4=mean(a$churisk_age[which(a$chbelief==-0.5&a$age1>25&a$yearstodate>5)])
a5=mean(a$churisk_age[which(a$chbelief==-1&a$age1>25&a$yearstodate>5)])
meano=c(a5,a4,a3,a2,a1)

b1=mean(a$churisk_age[which(a$chbelief==1&a$age1<26&a$yearstodate>5)])
b2=mean(a$churisk_age[which(a$chbelief==0.5&a$age1<26&a$yearstodate>5)])
b3=mean(a$churisk_age[which(a$chbelief==0&a$age1<26&a$yearstodate>5)])
b4=mean(a$churisk_age[which(a$chbelief==-0.5&a$age1<26&a$yearstodate>5)])
b5=mean(a$churisk_age[which(a$chbelief==-1&a$age1<26&a$yearstodate>5)])
meany=c(b5,b4,b3,b2,b1)


plot(meano,ylim=c(-0.5,1),type="l",pch=15,lwd=2,
     xlab="Change in Preference for Higher Social \n Spending since first wave",
     ylab="Mean change in unemployment \n risk since first wave",
     main="(c) Higher Social Spending and \n Unemployment Risk"
     ,cex.main=0.9,xaxt="n",cex.lab=0.8,las=2,cex.axis=0.9,font.main=1)
box()
axis(1,at=c(1,2,3,4,5),labels=c(-2,-1,0,1,2),cex.axis=0.9)
lines(meany,type="l",lty=2,lwd=2)
abline(h=0)


# Means, old
a1=mean(a$churisk_age[which(a$chbelieft==1&a$age1>25&a$yearstodate>5)])
a2=mean(a$churisk_age[which(a$chbelieft==0.5&a$age1>25&a$yearstodate>5)])
a3=mean(a$churisk_age[which(a$chbelieft==0&a$age1>25&a$yearstodate>5)])
a4=mean(a$churisk_age[which(a$chbelieft==-0.5&a$age1>25&a$yearstodate>5)])
a5=mean(a$churisk_age[which(a$chbelieft==-1&a$age1>25&a$yearstodate>5)])
meano=c(a5,a4,a3,a2,a1)

b1=mean(a$churisk_age[which(a$chbelieft==1&a$age1<26&a$yearstodate>5)])
b2=mean(a$churisk_age[which(a$chbelieft==0.5&a$age1<26&a$yearstodate>5)])
b3=mean(a$churisk_age[which(a$chbelieft==0&a$age1<26&a$yearstodate>5)])
b4=mean(a$churisk_age[which(a$chbelieft==-0.5&a$age1<26&a$yearstodate>5)])
b5=mean(a$churisk_age[which(a$chbelieft==-1&a$age1<26&a$yearstodate>5)])
meany=c(b5,b4,b3,b2,b1)

plot(meano,ylim=c(-0.5,1),type="l",lwd=2,
     xlab="Change in Preference for Higher Taxes \n on Rich since first wave",
     ylab="Mean change in unemployment \n risk since first wave",
     main="(d) Higher Taxes on the Rich and \n Unemployment Risk"
     ,cex.main=0.9,xaxt="n",cex.lab=0.8,las=2,cex.axis=0.9,font.main=1)
box()
axis(1,at=c(1,2,3,4,5),labels=c(-2,-1,0,1,2),cex.axis=0.9)
lines(meany,type="l",lty=2,lwd=2)
abline(h=0)

par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend(-0.45,-0.9, c("Respondents Aged 26+","Respondents Aged 18-25"), xpd = TRUE, horiz = TRUE,
       inset = c(0,0), lty=c(1,2),cex = 0.8)
dev.off()









